Grouping time series by pairwise measures of redundancy 
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A novel approach is proposed to group redundant time series in the frame of causality. It assumes 
, that (i) the dynamics of the system can be described using just a small number of characteristic 

^SJ ' modes, and that (ii) a pairwise measure of redundancy is sufficient to elicit the presence of correlated 

degrees of freedom. We show the application of the proposed approach on fMRI data from a resting 
human brain and gene expression profiles from HeLa cell culture. 
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• ^ , Over the last years the interaction structure of maiiy complex systems has been mapped in terms of graphs, which 
■ can be characterized using tools of statistical physics [l[ . Dynamical networks model physical and biological behavior 
in many applications; examples range from networks of neurons 0, Josephson junctions arrays Q to genetic networks 
protein interaction nets f^l and metabolic networks Q. Synchronization in dynamical networks is influenced by 
the topology of the network [3] . The inference of dynamical networks is related to the estimation, from data, of the 
^ flow of information between variables. Two major approaches are commonly used to estimate the information flow 

between variables, transfer entropy [sl] and Granger causality 
0^ ■ An important notion in information theory is the redundancy in a group of variables, formalized in jlQi] as a 
' generalization of the mutual information. A formalism to recognize redundant and synergetic variables in neuronal 
, ensembles has been proposed in 11 1 and generalized in [ip. Recently a quantitative definition to recognize redundancy 



\^ ' and synergy in the frame of causality has been provided jl3l | and it has been shown that the maximization of the total 
' causality, over all the possible partitions of variables, is connected to the detection of groups of redundant variables; the 
1 search over all the partitions is unfeasible but for small systems. We remark that the information theoretic treatments 
[ ' of groups of correlated degrees of freedom can reveal their functional roles in complex systems. The purpose of this 
^ work is to propose a simple approach to find groups of causally redundant variables (groups of variables sharing the 
same information about the future of the system), which can be applied also to large systems. The main assumption 
underlying our approach is that the essential features of the dynamics of the system under consideration are captured 
' using just a small number of characteristic modes. Hence we use principal components analysis to obtain a compressed 
representation of the future state of the system. Then, we introduce a pairwise measure of the redundancy w.r.t. 
the prediction of the next configuration of the modes, thus obtaining a weighted graph. Finally, by maximizing the 
modularity Q, we find the natural modules of this weighted graph and identify them with the groups of redundant 
variables. In the following section we describe the method. In section II we describe the application of the method to 
fMRI data, and in section III to a gene expression data-set. Some conclusions are drawn in section IV. 



I. METHOD 



Let us consider n time series {a;i(t)}i=i^..._„; after a linear transformation, we may assume all the time series to be 
normalized and with zero mean. The lagged times series are denoted Xi{t) = Xi{t — 1) . We make the hypothesis that 
the dynamics of the system under consideration may be described in terms of a few modes, and that these modes may 
be extracted by principal components analysis, as follows. Calling x the nxT matrix with elements Xi{t), we denote 
{'u-a{t)}a=i,...,nx the (normalized) eigenvectors of the matrix x^x corresponding to the largest nx eigenvalues. The 
T-dimensional vectors Uait) summarize the dynamics of the system; the lagged correlations of the system determine 
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to what extent the modes u may be predicted on the basis of the Xi{t) variables. 

PreUminarily, we select the variables which are significatively correlated with the modes u. For each i and each a 
we evaluate the probability pia that the correlation between Xi and Ua is due to chance, obtained by Student's t test. 
We compare pia with the 5% confidence level after Bonferroni correction (the threshold is 0.05/(n x n\)) and retain 
only those variables Xi which arc significatively correlated with at least one mode. The variables thus selected will 
be denoted {li(t)}i=i,...,Ar, N being their cardinality. 

The second step of the present approach is the introduction of a bivariate measure of redundancy, as follows. For 
each pair of variables Yi and Yj , we denote Pi the projector onto the one-dimensional space spanned by Yi and Pj the 
projector onto the space corresponding to Yj; Pij is the projector onto the bi-dimensional space spanned by Yi and 



Yj. Then, we define: 



J2{muc.\\^ + \\pjna\\^~mjuc.\\^); (1) 



according to the discussion in [13[, Cij is positive (negative) if variables Yi and Yj are redundant (synergetic) w.r.t. 
the prediction of the future of the system. In other words, if Yi and Yj share the same information about u, then Cij 
is positive. 

In the third step, the matrix Cij is used to construct a weighted graph of N nodes, the weight of each link measuring 
the degree of redundancy between the two variables connected by that link. By maximization of the modularity p3 | , 
the number of modules, as well as their content, is extracted from the weighted graph. Each module is recognized as 
a group of variables sharing the same information about the future of the system. 

As an example we report the following example. Let us consider the following autoregressive system: 

= 0.677t_i+ 0.1^1 

where ^ are i.i.d. unit variance Gaussian variables. By construction, if) is caused by 77 and viceversa. A system of 50 
time series is constructed as follows. For i — 1, . . . , 10: 

x^{t) =Vt+0.2pj, 

xw+,{t) = 7?t + 0.2pi°+\ 
e3 i n 0^20+1 



X2o+,(i) =e? + 0.2pf+\ (3) 
X3o+.(i) =et' + 0.2pf+\ 

l-L\ 40+i 



where p and ^ are i.i.d. unit variance Gaussian variables. Starting from a random initial configuration, the above 
equations are iterated and, after discarding the initial transient regime, rig consecutive samples of the system are 
stored for further analysis. Note that the first ten variables share the same information corresponding to ip, whilst 
the second ten variables share the information of rj. The variables Xi, with i — 21, . . . , 30, form a group of variables 
with correlations at equal times, similarly to the group of variables with i — 31,..., 40. The variables Xi, with 
i = 41, . . . , 50, correspond to pure noise. In figure ([T]) the equal-times correlations of the system are depicted, for a 
typical case with = 500, showing four groups of correlated variables. We perform the principal components analysis 
and retain a variable number nx of modes for the analysis. 

In figure ([2]), top-left, we depict N, the number of selected variables, as a function of n\. For nx = 3,4,5, twenty 
variables {xi with i = 1, . . . , 20) are selected; nineteen variables for nx = 1, 2, 6, 7, 8. Then, for each value of nx, the 
quantities Cij are evaluated. We find that, in this example, the matrix Cij is non- negative and can be treated as a 
weighted graph. 

In figure ([2]), top-right, we plot the number of modules Nm we find by applying the method described in 14!] to the 
matrix ; the method correctly recognizes the two modules for each value of Ux ■ 

In figure bottom-left, we plot a measure of the stability of the partition while going from — 1 to n^, defined 
as follows. We consider all the pairs of variables that are selected both in correspondence of tia — 1 and nx- The 
stability is one minus the fraction of pairs such that the variables are recognized to be in the same module in one 
instance and in different modules in the other instance. In this case the stability is always one; when the method is 
applied to real data, the stability curve may be helpful to fix the optimal number of modes nx- 

Finally, in figure (HJ, bottom-right, the eigenvalues of the matrix x^ x are depicted. In this case it is clear that the 
optimal number of modes is four. 

We remark that a suitable number of samples is needed to obtain reliable results. In figure ([3]) we depict the number 
of selected variables, for this example, as a function of for three choices of nx- it vanishes as the number of samples 
decreases. 
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II. MODULAR ORGANIZATION OF BRAIN ACTIVITY 

The £MRI signal can be regarded as a proxy for the underlying neural activity. Remote regions of the brain do not 
operate in isolation and there is a growing interest in studying the interactions and connectivity patterns between 
these regions, which have been investigated by independent component analysis [isj . principal components analysis 
and other approaches. Temporal and spatial functional networks, corresponding to spontaneous brain activity 
in humans, were derived in [l7| on the basis of the equal-time correlation matrix. Modularity in the resting state 
of the human brain has also been studied in [l8l - l20| . The connectivity structure of brain networks extracted from 
spontaneous activity signals of healthy subjects and epileptic patients has been analyzed in [2l|, . 

Here we consider fMRI data from a subject in resting conditions, with sampling frequency 1 Hz, and number of 
samples equal to 500. A prior brain atlas is utilized to parcellate the brain into ninety cortical and subcortical regions, 
and a single time series is associated to each region. All the ninety time series are then band-passed in the range 

0. 01-0.08 Hz. 

In figure (jj]), top-left, we depict N, the number of selected variables, as a function of n\. For n\ > 3, all the ninety 
regions are recognized as influencing the future of the system. For each value of tt-a, the quantities Cij are evaluated. 
In figure (|3]), top- right, the number of modules we find by applying the method described in to the matrix 
Cij is depicted; this plot suggests the presence of four modules for 4 < nx < 8. These values are the most stable, as 
it is clear from figure ([2]), bottom- left, where we plot the measure of the stability of the partition while going from 
n\ — 1 to n\. It follows that the optimal value of n\ is four, corresponding to a graph structure with four modules 
and modularity equal to 0.3. We find that module 1 includes brain regions from ventral medial frontal cortices which 
are primarily specialized for anterior default mode network, module 2 is typical referred to as posterior default mode 
network, module 3 mainly corresponds to executive control network and module 4 refers to the subcortical network. 
In figure bottom-right, the eigenvalues of the matrix x are depicted. 

It is worth comparing the histogram of the values of Cij, in this example (figure ([5]) -bottom), with those correspond- 
ing to a random choice of the modes u (figure ©-top). In the random case, the pairs of variables are either redundant 
or synergetic, and the typical values of c are very small. On the data set at hand, the magnitude of the values of c's 
is much higher and variables are mostly redundant; indeed the c's are negative (with small absolute value) only for a 
few pairs of regions. We remark that the presence of a few small and negative weights does not influence significantly 
the output from the modularity algorithm of fT^]: the output does not change if all c's with absolute value less than a 
threshold are set to zero (the threshold being chosen so that all the elements of the resulting matrix are nonnegative). 

Averaging the time series belonging to each module, we obtain four time series and we evaluate the causalities 
between them: the result is displayed in figure ©. It is interesting to observe that module 4 infiuences all the three 
other modules but is not influenced by them, it is an out-degree hub; this is consistent with the fact that it corresponds 
to subcortical brain. Another striking feature is the clear interdependencies between modules 2 and 3. The reliability 
of this pattern needs to be assessed on a large population of subjects. 

III. HELA GENE EXPRESSION REGULATORY NETWORK 

HeLa [1^ is a famous cell culture, isolated from a human uterine cervical carcinoma in 1951. HeLa cells have 
somehow acquired cellular immortality, in that the normal mechanisms of programmed cell death after a certain 
number of divisions have somehow been switched off. We consider the HeLa cell gene expression data of [l^l . Data 
corresponds to 94 genes and 48 time points, with an hour interval separating two successive readings (the HeLa 
cell cycle lasts 16 hours). The 94 genes were selected, from the full data set described in [2^, on the basis of the 
association with cell cycle regulation and tumor development. This data has been also considered in [2^. The static 
correlation analysis between time series, which is the result of regulation mechanisms with time scales faster than the 
sampling rate, revealed a highly related network with the presence of two modules: the flrst module was recognized as 
corresponding to the regulatory network of the transcriptional factor NFkB [27], whilst the second module appeared 
to be orchestrated by transcriptional factors p53 and STAT3. Use of bivariate Granger causality, in [2^, has put 
in evidence 19 causality relationships acting on the time scale of one hour, all involving genes playing some role in 
processes related to tumor development. 

As stated in [l^, fundamental patterns underlie gene expression proflles. This suggest the use of the proposed 
approach on gene expression time series. In flgure (jG)) we describe the application of the proposed approach on 
the HeLa data-set. The stable partition corresponds to ti\ = 4 and consists of two modules of 9 and 7 genes (the 
modularity is 0.1). The first module is characterized by the transcriptional factor NFkB and consists of NFkB, MCP- 

1, ICAM-1, Bcl-XL, lAP, A20, c-myc, TSPl, and Mcl-l. The second module is related to the transcriptional factor 
JunB, known to be a regulator of life and death of cells '29'], and consists of JunB, IL-6, IkappaBa, P21, Noxa, c-jun 
and NRBP. Averaging the time series belonging to each module, and evaluating the causality between the two time 
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series thus obtained, we obtain a relevant (0.145) causality of the first module on the second one. We note that all the 
sixteen genes selected by our approach were recognized as interacting in (2^ ; nine of them were involved also in the 
causalities described in [2^. It is not surprising that different methods, on the same data set, provide slightly different 
results: currently available data size and data quality make the reconstruction of gene regulatory networks from gene 
expression data a challenge. In figure ([SJ-bottom) we depict the histogram of the values of c on this data-set, with 
those corresponding to choosing randomly the modes u (figure ©-top). On this data, the values of c that we obtain 
are significatively greater than those one finds in the random case, and all the pairs are redundant. 

IV. CONCLUSIONS 

Grouping redundant time series reveals their functional roles in complex systems. In the frame of causal approaches, 
grouping redundant time series may reflect directed influence of one group over another. In this work we have proposed 
a novel approach which assumes that (i) the dynamics of the system can be described using just a small number of 
characteristic modes, and that (ii) a pairwise measure of redundancy is sufficient to elicit the presence of correlated 
degrees of freedom. Grouping is provided by the identification of the modules of the weighted graph of redundancies. 
The method may be seen as an alternative to the analysis of [l^, which can be applied also for large systems. We 
have shown the effectiveness of the proposed approach in two applications, the analysis of fMRI data and the analysis 
of gene expression data. In both these applications usually linear interactions are sought for, and only lags of 1 are 
considered, therefore here we limited to consider a linear pairwise measure of redundancy, and considered only lags 
of 1. The generalization of the pairwise redundancy measure, here introduced, to the nonlinear case and to lags of 
higher order is matter for further work, along the lines described in 13]. 
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FIG. 1: The correlation matrix of the simulated example, showing four groups of variables correlated at equal times. 



25 



20 



15 



2 4 6 8 10 



m 0.5 



2 4 6 8 10 



2 4 6 8 10 



10 20 30 



FIG. 2: (Top-left) Concerning the simulated example, the number of selected variables A'' is plotted versus nx, the number of 
modes. (Top-right) The number of modules, obtained by modularity maximization, of the matrix dj , whose elements measure 
the pairwise redundancy. (Bottom- left) The measure of the stability of the partition, going from nx — 1 to nx, is plotted versus 
nx- (Bottom-right) The eigenvalues of the matrix x^x are depicted. 
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FIG. 3: The number of selected variables, for the simulated example, is depicted as a function of the number of samples Us for 
riA = 2,3,4. 




FIG. 4: (Top-left) Concerning the fMRI application, the number of selected regions A'^ is plotted versus nx, the number of 
modes. (Top-right) The number of modules, obtained by modularity maximization, of the matrix aj , whose elements measure 
the pairwise redundancy. (Bottom- loft) The measure of the stability of the partition, going from nx — 1 to nx, is plotted versus 
nx- (Bottom-right) The eigenvalues of the matrix x'^ x are depicted. 
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FIG. 6: The causalities between the four modules of the fMRI application. 
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FIG. 7: (Top- left) Concerning the genetic application, the number of selected regions A'^ is plotted versus nx, the number of 
modes. (Top-right) The number of modules, obtained by modularity maximization, of the matrix dj, whose elements measure 
the pairwise redundancy. (Bottom-left) The measure of the stability of the partition, going from — 1 to nx, is plotted versus 
nx- (Bottom-right) The eigenvalues of the matrix x^a; are depicted. 




FIG. 8: The histogram of the values of the pairwise redundancy dj, in genetic application (bottom), and choosing randomly 
the modes u (top) 



